%load_ext cythonmagic
import numpy
import math
import scipy
import scipy.io
import scipy.misc
import scipy.cluster
import scipy.cluster.vq
import matplotlib
import matplotlib.pyplot
import os
import IPython
import IPython.parallel
import itertools
import random
import sklearn
import sklearn.decomposition
import sklearn.manifold
import sklearn.cluster
import sklearn.feature_extraction
import sklearn.feature_extraction.text
base_path = "/u/mlrobert/code/local/2013_drawing_assistant/data/"
img_path = base_path+"renderings/bunny2/apparent_ridges.num_latitude_lines=20.num_longitude_lines=20/"
filter_path = base_path+"filters/gabor/num_thetas=08/"
local_feature_cluster_centroids_path = \
base_path+"local_feature_cluster_centroids/bunny2/"+ \
"apparent_ridges.num_latitude_lines=20.num_longitude_lines=20.gabor.num_thetas=08.galif.patch_width=15.num_samples=32.num_tiles=04.num_samples=1000000.k=5000/"
global_features_path = \
base_path+"global_features/bunny2/"+ \
"apparent_ridges.num_latitude_lines=20.num_longitude_lines=20.gabor.num_thetas=08.galif.patch_width=15.num_samples=32.num_tiles=04.num_samples=1000000.k=5000/"
img_names_exts = sorted(os.listdir(img_path))
local_feature_cluster_centroids_path_name_ext = local_feature_cluster_centroids_path+"local_feature_cluster_centroids.mat"
global_features_names_exts = sorted(os.listdir(global_features_path))
num_thetas = len(sorted(os.listdir(filter_path)))
num_images = len(img_names_exts)
codebook_size = 5000
filter_path = base_path+"filters/gabor/num_thetas=08/"
tmp = []
for theta_index in range(num_thetas):
current_g_mat = scipy.io.loadmat(filter_path+"theta=%02d.mat" % theta_index)
current_g = current_g_mat["g"]
tmp.append(current_g)
g = numpy.array(tmp)
local_feature_cluster_centroids_mat = scipy.io.loadmat(local_feature_cluster_centroids_path_name_ext)
local_feature_cluster_centroids = local_feature_cluster_centroids_mat["local_feature_cluster_centroids"]
all_global_features = []
for global_features_name_ext in global_features_names_exts:
global_features_mat = scipy.io.loadmat(global_features_path+global_features_name_ext)
global_features = global_features_mat["global_features"].ravel()
all_global_features.append(global_features)
all_global_features_2d = numpy.array(all_global_features)
tf_idf_transformer = sklearn.feature_extraction.text.TfidfTransformer()
tf_idf_all_global_features_2d = numpy.array(tf_idf_transformer.fit_transform(all_global_features_2d).todense()).squeeze()
img_index = 124
img = scipy.misc.imread(img_path+img_names_exts[img_index]).astype(numpy.float32)[:,:,0]
img = (255.0 - img) / 255.0
figsize(4,4)
imshow(img, cmap="gray");
title(img_names_exts[img_index])
colorbar();
num_samples_width = 32
patch_width_as_fraction_of_width = 0.15
use_squared_filter_response = True
num_samples_height = int((float(img.shape[0])*float(num_samples_width))/float(img.shape[1]))
samples_y = numpy.linspace(0, img.shape[0]-1, num_samples_height+2).astype(numpy.int32)
samples_x = numpy.linspace(0, img.shape[1]-1, num_samples_width+2).astype(numpy.int32)
patch_size = int(img.shape[1] * patch_width_as_fraction_of_width)
if patch_size % 2 == 0:
patch_size = patch_size + 1
patch_half_size = int((patch_size - 1) / 2)
filter_responses_padded = zeros((num_thetas, img.shape[0]+patch_size-1, img.shape[1]+patch_size-1), dtype=float32)
def compute_filter_responses_padded():
for theta_index in range(num_thetas):
F_img = numpy.fft.fft2(img)
F_img_centered = numpy.fft.fftshift(F_img)
F_img_centered_filtered = F_img_centered * g[theta_index]
F_img_centered_filtered_uncentered = numpy.fft.ifftshift(F_img_centered_filtered)
filter_response = numpy.fft.ifft2(F_img_centered_filtered_uncentered)
if use_squared_filter_response:
filter_response = abs(filter_response)*abs(filter_response)
else:
filter_response = abs(filter_response)
filter_responses_padded[theta_index,patch_half_size:img.shape[0]+patch_half_size,patch_half_size:img.shape[1]+patch_half_size] = \
filter_response
compute_filter_responses_padded()
%%cython
import numpy
cimport numpy
import scipy
import scipy.cluster
import scipy.cluster.vq
ctypedef numpy.float32_t FLOAT32_DTYPE_t
ctypedef numpy.int32_t INT32_DTYPE_t
cdef inline float _mean(numpy.ndarray[FLOAT32_DTYPE_t, ndim=2] a):
cdef int y = 0
cdef int x = 0
cdef float s = 0
cdef float n = a.shape[0]*a.shape[1]
for y in range(a.shape[0]):
for x in range(a.shape[1]):
s = s+a[y,x]
return s / n
def compute_local_features(
numpy.ndarray[FLOAT32_DTYPE_t, ndim=2] img,
numpy.ndarray[FLOAT32_DTYPE_t, ndim=3] filter_responses_padded,
int num_thetas,
int num_samples_height,
int num_samples_width,
int patch_size):
cdef int num_tiles_per_dimension_per_patch
cdef float local_feature_norm_accept_threshold
num_tiles_per_dimension_per_patch = 4
local_feature_norm_accept_threshold = 0.01
cdef numpy.ndarray[INT32_DTYPE_t, ndim=1] samples_y
cdef numpy.ndarray[INT32_DTYPE_t, ndim=1] samples_x
cdef numpy.ndarray[INT32_DTYPE_t, ndim=1] tile_bounds
cdef numpy.ndarray[FLOAT32_DTYPE_t, ndim=2] patch
cdef numpy.ndarray[FLOAT32_DTYPE_t, ndim=2] tile
cdef numpy.ndarray[FLOAT32_DTYPE_t, ndim=1] tile_1d
cdef numpy.ndarray[FLOAT32_DTYPE_t, ndim=1] local_features_1d
cdef int sample_y, sample_x, theta_index, tile_index_y, tile_index_x, patch_index_y, patch_index_x
samples_y = numpy.linspace(0, img.shape[0]-1, num_samples_height+2).astype(numpy.int32)
samples_x = numpy.linspace(0, img.shape[1]-1, num_samples_width+2).astype(numpy.int32)
tile_bounds = numpy.linspace(0, patch_size, num_tiles_per_dimension_per_patch+1).astype(numpy.int32)
all_local_features_list = []
for sample_y in samples_y:
for sample_x in samples_x:
local_features_list = []
for theta_index in range(num_thetas):
filter_response_padded = filter_responses_padded[theta_index,:,:]
patch = filter_response_padded[sample_y:sample_y+patch_size, sample_x:sample_x+patch_size]
tile = numpy.zeros((num_tiles_per_dimension_per_patch, num_tiles_per_dimension_per_patch), dtype=numpy.float32)
for tile_index_y in range(num_tiles_per_dimension_per_patch):
for tile_index_x in range(num_tiles_per_dimension_per_patch):
tile[tile_index_y, tile_index_x] = _mean(patch[tile_bounds[tile_index_y]:tile_bounds[tile_index_y+1], \
tile_bounds[tile_index_x]:tile_bounds[tile_index_x+1]])
tile_1d = tile.ravel()
local_features_list.append(tile_1d)
local_features_1d = numpy.array(local_features_list).ravel()
if numpy.linalg.norm(local_features_1d) > local_feature_norm_accept_threshold:
all_local_features_list.append(local_features_1d)
return numpy.array(all_local_features_list, dtype=numpy.float32)
local_features = compute_local_features(img, filter_responses_padded, num_thetas, num_samples_height, num_samples_width, patch_size)
def compute_global_features():
local_features_codebook_indices, distances = scipy.cluster.vq.vq(local_features, local_feature_cluster_centroids)
hist, bin_edges = numpy.histogram(local_features_codebook_indices, bins=codebook_size)
tf_idf_hist = numpy.array(tf_idf_transformer.transform(hist).todense()).squeeze()
compute_global_features()
plot(tf_idf_hist);
def compute_distances():
cos_distances = dot(tf_idf_all_global_features_2d, tf_idf_hist)
distances = arccos(cos_distances)
sorted_indices = argsort(distances)
print "%s (distance = %f) " % (img_names_exts[sorted_indices[0]], distances[sorted_indices[0]])
print "%s (distance = %f) " % (img_names_exts[sorted_indices[1]], distances[sorted_indices[1]])
print "%s (distance = %f) " % (img_names_exts[sorted_indices[2]], distances[sorted_indices[2]])
print "%s (distance = %f) " % (img_names_exts[sorted_indices[3]], distances[sorted_indices[3]])
print "%s (distance = %f) " % (img_names_exts[sorted_indices[4]], distances[sorted_indices[4]])
plot(distances);
latitude=07.longitude=04.png (distance = 0.000000) latitude=08.longitude=04.png (distance = 0.405799) latitude=09.longitude=04.png (distance = 0.451262) latitude=07.longitude=05.png (distance = 0.583955) latitude=09.longitude=05.png (distance = 0.627951)
def find_matches():
compute_filter_responses_padded()
local_features = compute_local_features(img, filter_responses_padded, num_thetas, num_samples_height, num_samples_width, patch_size)
compute_global_features()
compute_distances()
%timeit -n2 -r3 find_matches()
2 loops, best of 3: 2.15 s per loop